#This is a replication file for survey experiment carried out in Sweden during July-August 2013.
#Kris-Stella Trump kstrump@gmail.com


library(foreign)

#Read in Uppsala data
UUdata <- read.csv("Exp3 - clean data UU.csv", na.strings="")
#Read in Kaninen data
SKdata <- read.csv("Exp3 - clean data kaninen.csv", na.strings="")

names(UUdata)
names(SKdata)

#Standardizing the datasets for a merge.
#The Uppsala sample did not have a question on status as student - add a column to reflect this
UUdata$studerande <- c(NA)
#Also add dummies that can identify the sample later if need be
UUdata$UU <- c(1)
UUdata$SK <- c(0)
SKdata$UU <- c(0)
SKdata$SK <- c(1)

#combine the datasets
alldata <- rbind(UUdata,SKdata)

#remove incompletes (participants who did not complete survey)
data <- alldata[alldata$V10==1,] 


#Create the variables needed for analysis

#Perceptions

data$percmin <- apply(data[,c("est_1_1_TEXT", "est_8_1_TEXT", "est_9_1_TEXT", "est_10_1_TEXT", "est_11_1_TEXT", "est_12_1_TEXT")], 1, min)
data$percmax <- apply(data[,c("est_1_1_TEXT", "est_8_1_TEXT", "est_9_1_TEXT", "est_10_1_TEXT", "est_11_1_TEXT", "est_12_1_TEXT")], 1, max)
data$percratio <- (data$percmax/data$percmin)
data$percindex <- log(data$percratio)
summary(data$percratio)


#Preferences total
data$prefmin <- apply(data[,c("tot.pref_x1", "tot.pref_x2", "tot.pref_x3", "tot.pref_x4", "tot.pref_x5", "tot.pref_x6")], 1, min)
data$prefmax <- apply(data[,c("tot.pref_x1", "tot.pref_x2", "tot.pref_x3", "tot.pref_x4", "tot.pref_x5", "tot.pref_x6")], 1, max)
data$prefratio <- (data$prefmax/data$prefmin)
data$prefindex <- log(data$prefratio)
summary(data$prefratio)

#Check for outliers
plot(data$prefratio)
plot(data$prefindex)
#Nothing extraordinary

#BJW indicator
data$bjw <- (data$bjw_1 + data$bjw_2 + data$bjw_3 + data$bjw_4 + data$bjw_5 + data$bjw_6 + data$bjw_7 + data$bjw_8) / 8
summary(data$bjw)
#Check distribution
plot(density(data$bjw, na.rm=T))

#Treatment indicator
#Information
data$info <- (is.na(data$tr.table))==FALSE
summary(data$info)

#Political affiliation
data$left <- data$party==1 | data$party==2 | data$party==11 
data$right <- data$party>11


#Descriptive statistics

#perceptions
summary(data$percratio)
sum(data$percratio>48) #6 people correct or higher 
6/250 #2.4% 

#preferences
summary(data$prefratio)

#bjw
summary(data$bjw)

#political affiliation
summary(data$left) #104/250 - 41%
summary(data$right) #67/250 - 27%

#age
hist(data$born)
summary(data$born) 
#mean 26, median 24, min 19, max 59

sum(data$gender==1) #150 female

sum(data$educ>=4, na.rm=T)
118/250 #=0.47
data$degree <- data$educ>=4



##### Balance table
balance <- zelig(info ~ born + gender + right + degree + bjw + percindex, model="ls", data=data)
summary(balance)

#Age (born+18 for age)
mean(data[data$info==0,"born"], na.rm=T)
mean(data[data$info==1,"born"], na.rm=T) 
#Gender
sum(data[data$info==0 & data$gender==1,"gender"], na.rm=T)/sum(data$info==0)
sum(data[data$info==1 & data$gender==1,"gender"], na.rm=T)/sum(data$info==1)
#BJW
mean(data[data$info==0,"bjw"], na.rm=T)
mean(data[data$info==1,"bjw"], na.rm=T) 
#Right wing
sum(data[data$info==0 & data$right==1,"right"], na.rm=T)/sum(data$info==0)
sum(data[data$info==1 & data$right==1,"right"], na.rm=T)/sum(data$info==1)
#College degree
sum(data[data$info==0 & data$degree==1,"degree"], na.rm=T)/sum(data$info==0)
sum(data[data$info==1 & data$degree==1,"degree"], na.rm=T)/sum(data$info==1)
#Perception of inequality
mean(data[data$info==0,"percindex"], na.rm=T)
mean(data[data$info==1,"percindex"], na.rm=T)



#Experimental results, Predicted values

#install.packages("Zelig")
library(Zelig)

#Show that treatment impacts preferences
z.out1 <- zelig(prefindex ~ info, model="ls", data=data)
summary(z.out1)


#Controlling for other relevant variables
z.out2 <- zelig(prefindex ~ info + bjw + percindex + right, model="ls", data=data)
summary(z.out2)


#Check impact of info on estimates
z.out.est <- zelig(percindex ~ info + bjw, model="ls", data=data)
summary(z.out.est)
#none -> successful randomization



################Graphing the results

library(Hmisc)

###Use model that only controls for experimental condition

#Create expected values. 
x.out <- setx(z.out1, info=c(0,1))
s.out <- sim(z.out1, x=x.out)
summary(s.out)

#Graph it

#order: 
#no info, info

mean <- c(1.1, 1.28) 
lci <- c(0.97, 1.15)
uci <- c(1.22, 1.41)
at <- c(3,7)

bjwresultslog <- as.data.frame(cbind(mean, uci, lci))
bjwresults <- exp(bjwresultslog)
bjwresults <- cbind(bjwresults,at)
mp <- barplot(bjwresults$mean,ann=F, ylim=c(1,10), xaxt="n", space=c(0.5),col=c("lightgrey","slategrey"))
bjwresults

pdf("inforesults_sweden.pdf")
barplot(bjwresults$mean,ann=F, ylim=c(1,6), xaxt="n", space=c(0.5),col=c("grey40", "lightgrey"),xpd=F, axes=F)
axis(1,at=c(0,3),lwd.ticks=0, labels=c("",""))
axis(2,at=c(1,10),lwd.ticks=1, labels=c("",""))
axis(2,at=c(2,3,4,5,6),lwd.ticks=1, labels=c(2,3,4,5,6))
errbar(mp,bjwresults$mean,bjwresults$uci,bjwresults$lci,add=T, pch=NA)
title(main="Recommended income ratios", ylab="Income ratio")
mtext(c("No information", "Information condition"), side=1, at=c(1,2.5), line=0)
dev.off()



#Should government reduce income differences? 
z.out.red1 <- zelig(opinion_7 ~ info, model="ls", data=data)
summary(z.out.red1)

z.out.red2 <- zelig(opinion_7 ~ info + bjw + percindex + right, model="ls", data=data)
summary(z.out.red2)


#Income differences are too large
z.out.diff <- zelig(opinion_6 ~ info, model="ls", data=data)
summary(z.out.diff)
#no impact

z.out.diff2 <- zelig(opinion_6 ~ info + bjw + percindex + right, model="ls", data=data)
summary(z.out.diff2)


#Trust politicians
z.out.trust <- zelig(trust.gov ~ info, model="ls", data=data)
summary(z.out.trust)
#no impact


